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Monte Carlo simulations within the grand canonical ensemble are used to explore the liquid- vapour 
coexistence curve and critical point properties of the Lennard-Jones fluid. Attention is focused on 
the joint distribution of density and energy fluctuations at coexistence. In the vicinity of the critical 
point, this distribution is analysed using mixed-field finite-size scaling techniques aided by histogram 
reweighting methods. The analysis yields highly accurate estimates of the critical point parameters, 
as well as exposing the size and character of corrections to scaling. In the sub-critical coexistence 
region the density distribution is obtained by combining multicanonical simulations with histogram 
reweighting techniques. It is demonstrated that this procedure permits an efficient and accurate 
mapping of the coexistence curve, even deep within the two phase region. 



I. INTRODUCTION 



The Lennard-Jones (LJ) fluid constitutes the prototype model for realistic atomic fluids and has been the focus of 
numerous simulation studies spanning well over 25 years [^|-pd|. The motivation for the longstanding interest in the 
model is its utility as a test-bed for ever more accurate and sophisticated theories of the liquid state. Contemporary 
theories ]T^-[T5|] now provide good agreement with simulation results over a wide range of noncritical temperatures. 
The continuing challenge, however, is to realise a similar degree of accuracy in the critical region, where the unbounded 
growth of correlations poses potentially serious difficulties for theory and simulation alike. 

One popular simulation method for studying the coexistence regime of fluid systems is the Gibbs ensemble Monte 
Carlo simulation technique (GEMC) introduced by Panagiotopoulos 0. In the GEMC method, the two coexisting 
phases separate into two physically detached, but thermodynamically connected boxes, the volumes of which are 
allowed to fluctuate under a constant pressure environment. Measurements of the particle density in each box provide 
estimates of the coexistence densities. Common practice is to fit the temperature dependence of these densities using 
a power law, the extrapolation of which yields estimates of the critical point parameters. The strength of the GEMC 
method lies in its elimination of the physical interface between the coexistence phases, the large free energy of which, 
plagues conventional grand canonical simulations of phase coexistence in the form of long lived metastable states and 
extended tunneling times. A large number of recent GEMC studies |16 testify to the method's efficacy in determining 
the subcritical coexistence properties of fluids. 

In the neighbourhood of the critical point, however, and for reasons discussed in references ]l7|-p0[|, the GEMC 
method cannot be relied upon to provide accurate estimates of the coexistence curve parameters. Instead it is 
necessary to employ finite-size scaling (FSS) techniques to probe the critical limit. FSS techniques were originally 
developed in the context of computer simulation studies of critical phenomena in spin models, and provide a highly 
effective route to infinite volume critical parameters from simulations of finite size [ pl|j22| . Recently their use has been 
extended to fluids by explicitly incorporating the consequences of the lack of symmetry between the coexisting phases 
]l7| . p3 25 1 . This reduced symmetry of fluids with respect to magnetic systems such as the Ising model, is manifest 
in the so-called 'field-mixing' phenomenon which is a crucial issue in the critical behaviour of fluids. The mixed-field 
FSS theory has been successfully employed in conjunction with simulations in the grand canonical ensemble (GCE) 
to study the critical behaviour of a number of critical fluid systems 17 2^ p6[ , including the two dimensional LJ fluid 
and a three dimensional lattice model for polymer mixtures. In the present work we extend these studies to the 3D 
LJ fluid making additional use of two important new methodological advances in computer simulation. 

The essential ideas underpinning the FSS methods we shall use, are the dual concepts of scale invariance and univer- 
sality. Precisely at criticality, the fluctuation spectra (distribution functions) of certain readily accessible observables 
assume scale invariant forms [^7|-^9|. Moreover these critical scaling functions are universal, being identical for all 
members of the same universality class. Experimental pp|-|32]|, theoretical |33|,[34| and simulation |2j| results show that 
the critical behaviour of simple fluids corresponds to the Ising universality class (to which all systems with short range 
interactions and a scalar order parameter belong). Scaling functions measured for the critical Ising model therefore 
constitute a hallmark of the Ising universality class, a fact that can be exploited to obtain accurate estimates of the 
critical point parameters of simple fluids. 

Besides the progress in extending the application of FSS concepts to fluid systems, two recent technical advances 
in computer simulation methods also greatly improve the efficiency with which one can tackle both the critical and 
subcritical regimes of model systems. The first is the histogram reweighting technique of Fcrrenberg and Swendsen 
pBj . This technique hinges on the observation that histograms of observables accumulated at one set of model 
parameters can be reweighted to yield estimates for histograms appropriate to another set of parameters. The 
histogram reweighting method has been found to be especially profitable close to the critical point where, owing to 
the large critical fluctuations, a single simulation affords reliable extrapolations over the entire critical region. 

The second technical advance came with the introduction by Berg and Neuhaus |36| of the multicanonical ensemble, 
use of which permits efficient Monte Carlo studies of two-phase coexistence even far below the critical point. The 
multicanonical method employs a prewcighting scheme to surmount the free energy barrier associated with formation 
of interfaces between the coexisting phases. This barrier grows rapidly as one moves further one from the critical 
point, quickly rendering conventional GCE simulations impractical. The multicanonical technique circumvents this 
problem by sampling not from a Boltzmann distribution, but from a preweighted distribution that is approximately 
flat between the coexisting densities. The desired Boltzmann distributed quantities are subsequently obtained by 
dividing out the preweighting factors from the measured histograms. Use of this method has been shown to reduce 
the tunnelling time to a simple power law in the system size fl37| . 

In the light of these new developments and the demand for ever increasing accuracy in estimates of the phase 
coexistence properties of prototype models, it seems appropriate to revisit the Lennard-Jones fluid with a view to 
performing a high precision study of its critical point and coexistence curve parameters. To this end we have carried 
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out a detailed simulation study of the model, bringing to bear all the aforementioned methodological advances. Our 
paper is organised as follows. We begin in section [II] by providing a short resume of the mixed-field FSS theory 
for the density and energy fluctuations of near-critical fluids. In section HI, we present measurements of the near 
critical scaling operator distributions as a function of system size. These distributions are analysed within the FSS 
framework, to yield extremely accurate estimates for the critical point and field mixing parameters. Turning then to 
the subcritical region we present the results of multicanonical simulations measurements of the coexistence density 
distributions. It is demonstrated how the measured coexistence densities can be used in conjunction with knowledge 
of the critical point parameters to construct the infinite volume coexistence curve. Finally, in section |iy| we detail 
our conclusions. 



II. THEORETICAL BACKGROUND 

In this section we provide a brief overview of the principal features of the mixed-field FSS theory of reference (jl7| . 

The system we consider is assumed to be contained in a volume L d (with d — 3 in the simulations to be considered 
below) and thermodynamically open so that the particle number can fluctuate. The observables on which we shall 
focus are the particle number density: 



L~ d N (2.1) 



and the dimensionless energy density: 



u = L- d {Aw)- 1 <S>{{r}) (2.2) 
where $({r}) is the configurational energy of the system which we assume takes the form: 

*({r})=5>(l r i- r jl)' (2.3) 

and where we assign the potential <j)(r) the familiar Lennard-Jones form 

4>(r) = 4w[{a/r) 12 - [a/rf] (2.4) 

with w the well-depth, and a a parameter that serves to set the length scale. 

Within the grand canonical ensemble (GCE), the joint distribution of density and energy fluctuations, pl{p, u), is 
controlled by the reduced chemical potential p and the well depth w (both in units of ksT). The critical point is 
located by critical values of the chemical potential \x c and well depth w c . Deviations of w and u from their critical 
values control the sizes of the two relevant scaling field that characterise the critical behaviour |58[ . In the absence of 
the special Ising ('particle-hole') symmetry, the relevant scaling fields comprise (asymptotically) linear combinations 
of the coupling and chemical potential differences |$9| : 

r = w c — w + s(/i — fi c ) h = fi — n c + r(w c — w) (2-5) 

where r is the thermal scaling field and h is the ordering scaling field. The parameters s and r are system-specific 
quantities controlling the degree of field mixing. In particular r is identifiable as the limiting critical gradient of the 
coexistence curve in the space of /i and w. The role of s is somewhat less tangible; it controls the degree to which 
the chemical potential features in the thermal scaling field, manifest in the widely observed critical singularity of the 
coexistence curve diameter of fluids Q . 

Conjugate to the two relevant scaling fields are scaling operators M. and £, which comprise linear combinations of 
the particle density and energy density p7| : 

M = — - — [p - su] £ = — - — [u - rp] (2.6) 

1 — sr 1 — sr 

The operator M. (which is conjugate to the ordering field h) is termed the ordering operator, while £ (conjugate to 
the thermal field) is termed the energy-like operator. In the special case of models of the Ising symmetry, (for which 
s = r = Q), M. is simply the magnetisation while £ is the energy density. 

The joint distribution of density and energy is simply related to the joint distribution of mixed operators: 

Pl(p,u) = -^— Pl (M 1 £) (2.7) 
1 — sr 
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Near criticality, and in the limit of large system size, Pl(M,£) is expected to be describable by a finite-size scaling 
relation of the form (ItJ : 



where 



and 



PL (M,£) s A+ A+^ i£ (A+ SM, A+5£, A M h, A £ t) 



A £ = asL 1 /", A M = a M L d -P/», A M A M = A £ A+ = L d 



5M = M- <M> C 



5£ = £- < £ > r 



(2.8a) 
(2.8b) 
(2.8c) 



The subscripts c in equations 2.8c signify that the averages are to be taken at criticality. Given appropriate choices for 
the non-universal scale factors cim and as (equa tion 2.8bj ), the function pm,s(A m SA4, A £ S£, Aj^h : A £ t) is expected 
to be universal. Precisely at criticality, equation |2.8a implies simply 



PL (M,£) ~ A M A+p*(A M 6M,A+S£) 



(2.9) 



where PM,s(x,y) = pM,s{x,y,0,0) is a function describing the universal and statistically scale invariant operator 
fluctuations characteristic of the critical point. 

For smaller system sizes, one anticipates that corrections to scaling associated with finite values of the irrelevant 
scaling fields will become significant | 3gj] . These irrelevant fields take the form a\T 8 + a^T 16 + ■ ■ ., where 9 is the 
universal correction to scaling exponent, whose value has been estimated to be 9 w 0.54 for the 3D Ising class [^o| . 
Incorporating the least irrelevant of these corrections into equation 2.E, one finally obtains 



p L {M,£)^A\ i A+p* M>£ {A + M 8M,Ap£,a 1 L- e l v ) 



(2.10) 



As we shall show, it is necessary to take account of such correction terms if highly accurate estimates of the critical 
parameters are to be obtained. 



III. MONTE CARLO STUDIES 



A. Computational details 

The Monte-Carlo simulations described here were performed using a Metropolis algorithm within the grand canoni- 
cal ensemble. The algorithm employed is similar in form to that described by Adams |4l],|^], but differs in the respect 
that only particle transfer (insertion and deletion) steps were implemented, leaving particle moves to be performed 
implicitly as a result of repeated transfers. Physically this choice is motivated by the need to direct the computational 
effort at the density fluctuations, which are the bottleneck for phase space evolution at coexistence. 

As is common practice in simulations of systems whose interparticle potential decays rapidly with particle separation, 
the Lennard- Jones potential was truncated in order to reduce the computational effort. In accordance with most 
previous studies of the LJ system, the cutoff radius was chosen to be r c — 2.5er, and the potential was left unshifted. 
It should be noted, however, that the choice of cutoff can have quite marked effects on the critical point parameters, 
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In order to facilitate efficient computation of interparticle interactions, the periodic simulation space of volume L 3 
was partitioned into m 3 cubic cells, each of side r c . This strategy ensures that interactions emanating from particles 
in a given cell extend at most to particles in the 26 neighbouring cells. We chose to study a range of system sizes 
corresponding to m = 3,4,5,6 and 7, containing at coexistence, average particle numbers of 135,320,625, 1080 and 
1715 respectively. For the m = 3,4 and 5 system sizes, equilibration periods of 10 5 Monte Carlo transfer attempts per 
cell (MCS) were utilised, while for the m = 6 and m — 7 system sizes up to 2 x 10 6 MCS were employed. Sampling 
frequencies ranged from 15 MCS for the m — 3 system to 50 MCS for the m = 7 system. The total length of the 
production runs was also dependent upon the system size. For the m = 3 system size, 1 x 10 7 MCS were employed, 
while for the m = 7 system, runs of up to 5 x 10 7 MCS were necessary. In the subcritical coexistence region, (studied 
using multicanonical simulations), runs of length 5 x 10 6 MCS were utilised. In both the sub-critical and critical 
coexistence regimes, the average acceptance rate for particle transfers was approximately 25%. 

In the course of the simulations, the observables recorded were the particle number density p and the energy density 
u. The joint distribution pl(p,u) was accumulated in the form of a histogram. In accordance with convention, we 
express p and u in reduced units: 
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per 



(3.1) 



We also n ote for future reference that the algorithm actually utilises not the true chemical potential p featuring in 
equation 2.5, but an effective chemical potential /i* to which the true chemical potential is related by 



At = n* + Mo - In (N/L d ) 



(3.2) 



where fio is the chemical potential in the non-interacting (ideal gas) limit. It is this effective value that features in 
the results that follow. 



B. The critical limit 



The most recent Gibbs ensemble simulation studies of the LJ fluid, (using r c = 2.5a), place the critical temperature 
at T* = 4/w = 1.176(8) |^3|. Using this estimate, we attempted to locate the liquid-vapour coexistence curve by 
performing a series of very short runs for the m = 4 system size, in which the effective chemical potential p.* was tuned 
until the density distribution exhibited a double peaked structure. Having obtained, in this manner, an approximate 
estimate of the coexistence chemical potential, a longer run comprising 2 x 10 7 MCS was performed to accumulate 
better statistics. Histogram reweighting was then applied to the resulting histogram enabling exploration of the 
coexistence curve in the neighbourhood of the simulation temperature. 

To facilitate a precise identification of the coexistence chemical potential, we adopted the criterion that the ordering 
operator distribution pl(A4) = J d£pL(M, £) must be symmetric in A4 — (M). This criterion is the counterpart of 
the coexistence symmetry condition for the Ising model magnetisation distribution. By simultaneously tuning ji and 
s in the reweighting of the joint distribution pl(p, u), estimates for the coexistence chemical potential and the value 
of the field mixing parameter s that satisfy this symmetry condition, were readily obtained. 

To obtain a preliminary estimate of the critical point parameters, the universal matching condition for the or- 
dering operator distribution Pl(M.) was invoked. As observed in sections [| and [il], fluid-magnet universality im- 
plies that the critical fluid ordering operator distribution pl(/A), must match the universal fixed point function 
Pm( x ) = / P*M e( x ^y)^y appropriate to the Ising universality class. The latter function is identifiable as the critical 
magnetisation distribution of the Ising model, the form of which is independently known from detailed simulation 
studies of large Ising lattices pij . Leaving aside for present the question of corrections to scaling, the apparent 
critical point of the fluid can thus be estimated by tuning the temperature, chemical potential and field mixing 
parameter s (within the reweighting scheme) such that pl(M) collapses onto p* M {x) . The result of applying this 
procedure for the m = 4 data set is displayed in figure [l] where the data has been expressed in terms of the scaling 
variable x = a]^L^I v {M - M c ) . The accord shown corresponds to a choice of the apparent critical parameters 
T*(L) = 1.1853(2), n%(L) = -2.7843(3). 

Using these estimates of the critical parameters, extensive simulations were then performed for each of the 5 system 
sizes m = 3 — 7 in order to facilitate a full finite-size scaling analysis. Reweighting was again applied to the resulting 
histograms to effect the matching of pl(-M) to p* M {x) thus yielding values of the apparent critical parameters. 
Interestingly, however, the apparent critical parameters determined in this manner were found to be L- dependent. 
The reason for this turns out to be significant contributions to the measured histograms from corrections to scaling, 
manifest as an L-dependent discrepancy between the critical operator distributions and their limiting Ising forms. In 
the case of the ordering operator distribution pl(M) , the symmetry of the Ising problem implies that the correction 
to scaling function is symmetric in M. — {M). In attempting to implement the matching to p* M (x) we therefore 
necessarily introduce an additional symmetric contribution to Pl(M) associated with a finite value of the scaling 
field t. This latter contribution has, coincidentally, a functional form that is very similar to that of the correction to 
scaling function, a result which of course make the cancellation of contributions possible. It follows, therefore, that 
the magnitude of two contributions must be approximately equal. 

Notwithstanding the added complications that corrections to scaling engender, it is nevertheless possible to extract 
accurate estimates of the infinite- volume critical parameters from the measured histograms. The key to accomplishing 
this is the known scaling behaviour of the corrections to scaling which, (recall equation 2.1C| ), die away with increasing 



system size like L~®l v . Now, since contributions to pl(-M) from finite values of r, grow with system size like \t\L x I v , 
it follows that implementation of the matching condition leads to a deviation of the apparent critical temperature 
T*{L) from the true critical temperature T* which behaves like 

T*(oo) - T*{L) oc L- [e+1)/v (3.3) 

In figure || we plot the apparent critical temperature T*(L) as a function of One observes that the data 

is indeed well described by a linear dependence, the least squares extrapolation of which yields the infinite-volume 
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estimate T* = 1.1876(3). The associated estimate for the critical chemical potential is fi* — —2.778(2). We note 
however, that although the coexistence value of /i* is tightly tied to T*, estimates of /i* are not directly affected by 
corrections to scaling in pl(M), since the function 5p l{M.)/ '8 (i is (to leading order) antisymmetric in M — (M) [ fl7J . 

Having acquired accurate estimates for the infinite volume values of T* and /i*, it is instructive to examine more 
closely the size and character of corrections to scaling in the operator distributions. Addressing first the ordering 
operator distribution, we show in figure §(a), the critical point form of pl{M) (expressed in terms of the scaling 
variable x = a^L^/ u {M — -M c )), for the two system sizes m — 4 and m = 7. Also shown is the universal fixed point 
function p* M (x) appropriate to the 3D Ising universality class. The corrections to scaling, manifest in the discrepancy 
between the fluid finite-size data and the limiting form, are clearly evident in the figure, especially for the m = 4 
system size. We note further that their form is qualitatively similar to those observed in the 2D Ising universality 
class @. 

A similar situation pertains to the energy operator distribution pl{£) — J dM.p l{M. , £■) ■ Figure ||(b) shows the 
form of this function, together with the limiting fixed point function p* £ (y) — J p* M £ {x,y)dx , identifiable as the 
critical energy distribution of the Ising model, and independently known from detailed Ising model studies |^5| . The 
data have all been expressed in terms of the scaling variable y = a^ 1 L d ~ l •• ' l '(£ — £ c ). One observes that in this case, 
the corrections to scaling are noticeably larger than for pl{M) , a fact that presumably reflects the relative weakness 
of critical fluctuations in £ compared to those in M . 

Turning now to the critical point field mixing parameters, s and r, the values of these quantities were assigned (as 
described in detail in reference pEf ) such as to optimise the mapping of the critical operator distributions onto their 
limiting fixed point forms, (cf. figure ||). The resulting estimates were, however, found to be slightly L-dependent for 
the smaller system sizes, an observation that may indicate a finite-size dependence of the scaling fields themselves [ fl6| . 
For the two largest system sizes this L-dependence is, however, small and we estimate s = —0.11(1), r = —1.02(1). 

The measured histograms also serve to furnish estimates of the exponent ratios (3/v and 1/v characterising the two 
relevant scaling fields h and r . These exponent ratios are accessible via the finit e-size scaling behaviour of pl(M) 
and Pl(£) at the critical point. Specifically, consideration of the scaling form |2.9| shows that the typical size of the 
critical fluctuations in the energy-like operator will vary with system size like S£ ~ i-Cd-V"), while the typical 
size of the fluctuations in the ordering operator vary like SA4 ~ L~^l v . Comparison of the standard deviation of 
these distributions as a function of system size thus affords estimates of the appropriate exponent ratios. In order 
to minimise systematic errors resulting from corrections to scaling, we have performed this comparison only for the 
two largest system sizes m — 6 and m — 7. From the measured variance of pl(M) for these two systems, wc find 
(3/v = 0.521(5), an estimate which compares very favourably with the three dimensional (3D) Ising estimate [[|7j of 
(3/v = 0.518(7). Given though that no allowances were made for corrections to scaling, the quality of this accord is 
perhaps slightly fortuitous. 

Carrying out an analogous procedure for pl{£) yields the estimate 1/v = 1.67(7), which does not agree to within 
error with the 3D Ising estimate 1/v = 1.5887(4). Here, though, we believe that the bulk of the discrepancy is 
traceable to the high sensitivity of Pl(£) with respect to the designation of the field mixing parameter r implicit in 
the definition of £ (cf. equation |2.6[ ). In the presence of sizeable corrections to scaling, it is somewhat difficult to 
gauge very accurately the infinite volume value of r from the mapping of Pl{£) onto p £ (y) ■ Studies of significantly 
larger system sizes than considered here would be necessary to alleviate this problem. 

Addressing now the critical density and energy distributions, figure ^ shows the measured forms of pl(p) and Pl(u) 
at the designated critical parameters. Clearly these distributions are to varying degrees asymmetric, a fact which 
(as explained in detail in reference p5| ) stems from field mixing effects. These field mixing contributions, (which are 
not be confused with corrections to scaling) die away with increasing L so that the limiting forms of both pl(p) and 
Pl(u) match the fixed point ordering operator distribution p* M {x). The approach to this limiting behaviour is indeed 
quite evident in figure [l| We note however, that the limiting form of the fluid critical energy distribution differs from 
that of the Ising model where limi^oo pl(u) = p* £ (y). This radical alteration to the limiting behaviour manifests the 
coupling that occurs in asymmetric systems between the ordering operator and energy-like operator fluctuations, the 
former of which dominate for large L ]25|,^6| . As a consequence one finds that for critical fluids the specific heat 

C v = L d ((u 2 ) - {u) 2 )/k B T 2 ~ V' v (3.4) 

in stark contrast to the behaviour in the Ising model for which C v ~ L"'". To recapture the Ising behavior it is 
instead necessary to consider the fluctations of the fluid energy-like operator £: 

L d ((£ 2 ) - (£) 2 )/k B T 2 - L a / V (3.5) 

As a further important consequence of field mixing, it transpires that measurements of the density and energy 
distributions at the infinite volume critical point do not afford direct estimates of the infinite volume critical density 
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and energy density. This was demonstrated in reference P5[ , where it was shown that the presence of field mixing 
contributions to pl(p) and pl(u) introduces a finite-size shift to their average values which behaves like the Ising 
energy: 

(p) c (L) - (p) c (co) ~ L-^- 1 '^ (3.6a) 
(u) c (L) - (u) c (oo) ~ L-^- 1 /^ (3.6b) 

Thus in order to obtain infinite volume estimates of p c and u c , it is necessary to perform a finite-size extrapolation of 
(p) c and (u) c to L — oo. In figure]^ we plot the values of (p) c and (it} C) corresponding to the distributions of figure |J, 
as a function of L - ^ -1 /^). Although no allowances have been made for corrections to scaling (the effects of which 
are certainly much smaller than those of field mixing), the data exhibit within the uncertainties, a rather clear linear 
dependence. Least-squares fits to the data yield the infinite volume estimates p* — 0.3197(4), u* — —0.187(2). 

We round off this subsection by summarising our results for the critical point parameters of the LJ fluid with 
r c = 2.5cr: 

T* = 1.1876(3), p* = -2.778(2) 
p* = 0.3197(4), u* = -0.187(2) 
s = -0.11(1), r = -1.02(1) (3.7) 

A comparison of these estimates with those of previous studies features in our concluding section. 



C. The subcritical coexistence region 

As described in section [j], conventional GCE simulation studies of the two-phase subcritical region encounter serious 
problems due to the large free energy barrier separating the coexisting phases. This can lead to pronounced metastabil- 
ity effects and protracted tunneling times between the phases. The multicanonical ensemble approach |3^] ameliorates 
these difficulties by artificially enhancing the frequency with which a simulation samples the interfacial configurations 
of intrinsically low probability. This enhancement is achieved by sampling not from a simple Boltzmann distribution 
with Hamiltonian H({r}, p), but from a modified distribution with effective Hamiltonian W({r}, p) — Tt({r}, p)+g(p), 
where g(p) is a preweighting function the specification of which is described below. For the case of the density, the 
preweighted distribution takes the form: 



N=L d p 



*<p)=t> n {/*■«} 



-l<S,({r}) +fJ ,L d p+g(p)] 



(3.8) 



where Z' is the multicanonical partition function, which is defined by equation 3.8. 

If one now chooses the preweighting function such that g(p) w lnp(p), where p(p) is the desired Boltzmann 
density distribution, one readily sees that p'(p) ~ constant \/p. To the extent that this condition is satisfied, the 
density thus performs a ID random walk over its entire domain, thereby allowing extremely efficient accumulation of 
the preweighted histogram p'(p). Once this histogram has been obtained, the desired Boltzmann weighted density 
distribution is regained as simply 

p(p)=p'(p)e-^ (3.9) 

Clearly for this approach to succeed, one requires a prior estimate of the function p{p) to use as the preweighting 
function. But p(p) is, of course, just the function we are trying to find! While feasible iterative schemes exist for 
estimating a suitable weight function f48fl , for the purposes of determining the coexistence curve distributions the task 
is considerably more straightforward. Knowledge of a near critical coexistence density distribution (easily obtainable 
since the free energy barrier to tunneling is small near T c ) can be used in conjunction with histogram reweighting 
and the equal peak- weight criterion |49]| to estimate the form of pl(p) for some other chosen point further down the 
coexistence line. The extrapolated estimate of the density distribution may then be used as the weight function in a 
multicanonical simulations at this new coexistence state point, yielding a new coexistence density distribution. The 
procedure is then simply repeated: histogram extrapolation of the new distribution being used to predict the weight 
function and coexistence parameters T*,p* for another state point still deeper into the subcritical region. In this 
manner one can systematically track along the coexistence curve in the space of p* and T* , obtaining at the same 
time the spectrum of coexistence density distributions. 
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We have implemented this strategy for the to = 4 system, employing the measured near-critical density distribution 
as our starting point. It was found that the histogram reweighting affords reliable extrapolations over rather a large 
temperature range: only 7 multicanonical simulations were required to reach temperature T* = 0.8T*. The resulting 
coexistence density distributions (corresponding to those temperatures at which the multicanonical simulations were 
actually performed) are depicted in figure ||. We note that for the lowest temperature studied, T* « 0.94, the ratio 
between the peak and trough heights of pl(p), is some 30 orders of magnitude! Such a difference would, of course, 
constitute an insurmountable barrier for a conventional grand canonical simulation. 

Away from the immediate vicinity of the critical point (where the correlation length £ <C L), the peak positions 
of the coexistence density distributions are expected to correspond to the densities of the infinite volume coexisting 
phases. This fact can be utilized in conjunction with the previously determined critical parameters to estimate the 
density-temperature phase diagram of the LJ fluid. In table Q we list the peak positions of our measured density 
distributions, which are also plotted as a function of T* in figure [?]. For comparison, the GEMC simulation data of 
Panagiotopoulos |l3) are also included. 

We have attempted to fit our non-critical density data to a power law of the form 



P± 



a\T* - T*\ ± b\T* - T*f (3.10) 



with T* and p* assigned the values given in equation |3.7| , and the order parameter exponent assigned the Ising 
estimate [3 — 0.3258 [^7j . The results of this fit are included in figure [7] (solid line) and correspond to a choice of 
the critical amplitudes a = 0.1824(3), b = 0.5226(4). As one observes, the data well away from the criti cal point are 



indeed very well fitted by the assumed form, implying both that the validity of the scaling form eq. |3.1C extends well 



into the subcritical region (a result also observed in many other simple fluids |lq|), and that the coexistence diameter 
singularity is undetectably small on the scale of our measurements. One further sees from figure ^ that for this system 
size (to = 4), systematic finite-size effects become apparent for T* > 0.95T C . Any power law fit to the density data 
that attempted to extrapolate to criticality by including datapoints closer to criticality than this, would thus run the 
risk of seriously overestimating the critical temperature Jl7]Jl8|,^5| . 

Finally, in this section we plot the coexistence curve in the space of p* and T* as obtained from the multicanonical 
simulations. Figure || shows this curve, together with the estimated critical point and the measured directions of the 
relevant scaling fields. 



IV. CONCLUSIONS 



In summary, we have employed recently developed mixed-field FSS techniques and histogram extrapolation methods 
to obtain highly precise estimates for the critical point parameters of the truncated and unshifted LJ fluid with 
r c = 2.5cr. Our measurements enable us to pinpoint the infinite volume critical temperature to within an uncertainty 
of 0.03%, considerably better than the accuracy of 1% (or more) typically quoted for other commonly used simulation 
techniques. 

Two recent studies have also reported values for the critical parameters of the LJ fluid with r c = 2.5a. That of Finn 
and Monson (§] corrected the equation of state data of Nicolas et al || for the discontinuity at r c and the absence of 
a long tail. Their resulting esti mat e of the critical temperature is T* — 1.23, which is clearly too high with respect 
to our own estimate (equation 3.7). Their value for the critical density p* = 0.32 does, on the other hand, agree 
well with our result, although since no error bars were quoted it is impossible to tell to what extent the accord is 
meaningful. 

By comparison, the GEMC simulation estimates of Panagiotopoulos T* = 1.176(8), p* — 0.33(1) correspond 
rather more closely to our results. In this GEMC study, a power law fit was made to subcritical coexistence density 
data, ignoring data points close to the critical point which are most influenced by finite-size effects. Although this 
approach certainly seems to reduce the systematic overestimate of T* that can occur in GEMC simulations when 
fitting all the available density data, it is not clear how much near-critical data should be discarded when the location 
of the critical point and the extent of finite-size effects are not known beforehand. Perhaps as a consequence of this, 
(as well as the neglect of corrections to scaling), the error bar on the value of T* quoted in jl3| does not overlap with 
ours. 

Turning now to the general computational issues raised by the present study, we have seen the great utility of FSS 
methods for probing the critical point region of fluids. The power of FSS techniques was also previously demonstrated 
in a related GCE study of the 2D LJ fluid jlT] . One significant drawback of this previous investigation, however, was 
its high computational cost. Histogram reweighting was not employed and consequently several long simulations were 
required for each L, in order to accurately locate the near-critical coexistence curve and the critical point. This in 
turn entailed the use of long runs on high performance parallel computers. By contrast, use of histogram reweighting 
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in the present work allowed the study of systems containing up to 5 times as many particles as those of the 2D study, 
while exercising the capabilities of only a pair of middle-range workstations. We thus believe that the combined use 
of FSS methods and histogram reweighting techniques, as espoused here, brings high precision studies of fluid critical 
phenomena within the reach of almost every pocket. 

The benefits offered by the use of multicanonical preweighting for simulations of the subcritical two-phase coexistence 
region are similarly impressive. In the previous GCE study of the 2D LJ fluid [0, phase coexistence could not be 
studied below about 0.98T C , due to the large free energy barrier separating the coexisting phases. However, by 
incorporating multicanonical preweighting into GCE simulations, we have seen that it is possible to probe much 
smaller subcritical temperatures with ease Q. 

Finally we remark that the techniques deployed here are not restricted to simple fluids, but can be combined with 
configurational bias Monte Carlo methods 51 5^] to facilitate accurate investigations of phase coexistence and critical 
phenomena in polymer systems. We also believe it should be feasible to apply the present method to systems with 
long-ranged interactions, e.g. coulombic fluids, although the computational workload would naturally increase rapidly 
with the range of the potential. We hope to report on such extensions in future work. 
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x=a M 'L p,v (M-M c ) 

FIG. 1. The measured form of the ordering operator distribution pl(M) for the m = 4 system size at the apparent critical 
parameters T* = 1.1853, fi* = —2.7843. Also shown for comparison is the universal fixed point ordering operator distribution 
Pm(x) . The data has been expressed in terms of the scaling variable x = aJ^L 13 ^ '(M—M c ), with the value of the non-universal 
scale factor aj^ chosen so that the distributions have unit variance. Statistical errors do not exceed the symbol sizes. 
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FIG. 2. The apparent reduced critical temperature, (as defined by the matching condition described in the text), plotted as 
a function of L~( e+1 ^ v , with 8 — 0.54 and v = 0.629 po| , |47| . The extrapolation of the least squares fit to infinite volume yields 
the estimate T* = 1.1876(3). 
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FIG. 3. (a) The ordering operator distribution p L (M) for the two system sizes m = 4 and m — 7 at the assigned critical 
parameters T* , fi*, expressed as function of the scaling variable x = aJ^L /3 ^ l/ (M — M c )- Also shown (solid line) is the universal 
fixed point ordering operator distribution p* M (x). (b) The energy operator distribution Pl(£) for the two system sizes m — 4 
and m = 7 at T c *,/ic ,expressed as a function of the scaling variable y = a^ 1 L d ~ 1 / 1 ' _ g c y Also shown (solid line) is the 
universal fixed point energy operator distribution p%{y) . In both cases the values of the non-universal scale factors a^ 1 or aT^ 
have been chosen to yield unit variance. 
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FIG. 4. (a) The density distribution at T* , fi* for the system sizes m = 4 — 7, (b) The corresponding energy density 
distributions. The lines are merely guides to the eye. Statistical errors do not exceed the symbol sizes. 
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FIG. 5. (a) The measured average density (p} c (L) at the designated critical point, expressed as a function of £-("*- V*). The 
least-squares fit yields an infinite volume estimate p c = 0.3197(4). (b) The measured average energy density (u) c (L) at the 
critical point, expressed as a function of L~ ( - d ~ 1 ^' 1 . The least-squares fit yields an infinite volume estimate u c = —0.187(2). 
In both cases we took 1/u = 1.5887 0. 
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FIG. 6. (a) Estimates of the coexistence density distributions for the m = 4 systems size, for a range of subcritical temper- 
atures, obtained as described in the text. The lines are merely guides to the eye. Statistical errors do not exceed the symbol 
sizes. 
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FIG. 7. The peak densities (filled triangles) corresponding to the distributions of figure q, plotted as a function of the reduced 
temperature. The coexistence diameter is also marked (filled squares). Statistical errors do not exceed the symbol sizes. Also 
shown (circles) are the Gibbs ensemble estimates of Panagiotopolous [M for a system of size L = 12a. The solid line represents 
a fit through T* , p* of the form p±- p c = a\T* -T*| ± b\T* - T*f, with a = 0.1824(3), b = 0.5226(4) and /3 = 0.3258 @. 
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FIG. 8. The line of liquid-vapour phase coexistence in the space of /x* and T*, for temperatures in the range 0.95 < T* < T* . 
The results were obtained by implementing the equal peak weight criterion for the density distribution [^9j in conjunction with 
the multicanonical simulations and histogram reweighting. Also shown are the measured directions of the relevant scaling fields. 
Statistical errors do not exceed the symbol sizes. 
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TABLE I. The peak densities corresponding to the coexis- 
tence curve distributions depicted in figure 



Temperature 


pv 


Pa 


1.1696 


0.1635(10) 


0.4855(10) 


1.1494 


0.1375(9) 


0.5155(9) 


1.1111 


0.1035(9) 


0.5599(9) 


1.0667 


0.0780(9) 


0.6031(9) 


1.0256 


0.0580(9) 


0.6380(9) 


0.9877 


0.0445(9) 


0.6665(9) 


0.9412 


0.0335(8) 


0.6915(8) 



